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ABSTRACT 

This paper demonstrates that rapid fracture of ideal brittle lattices nat- 
urally involves phenomena long seen in experiment, but which have been 
hard to understand from a continuum point of view. These idealized mod- 
els do not mimic realistic microstructure, but can be solved exactly and 
understood completely. First it is shown that constant velocity crack so- 
lutions do not exist at all for a range of velocities starting at zero and 
ranging up to about one quarter of the shear wave speed. Next it is shown 
that above this speed cracks are by and large linearly stable, but that at 
sufficiently high velocity they become unstable with respect to a nonlinear 
micro- cracking instability. The way this instability works itself out is re- 
lated to the scenario known as intermittency, and the basic time scale which 
governs it is the inverse of the amount of dissipation in the model. Finally, 
we compare the theoretical framework with some new experiments in Plex- 
iglas, and show that all qualitative features of the theory are mirrored in 
our experimental results. 

PACS: 62.20.Mk, 46.30.Nz Submitted to Journal of the Mechanics and Physics of Solids 



Introduction 



The classic theory of fracture (FREUND 1990; KANNINEN AND POPELAR 1985) treats 
cracks as mathematical branch cuts which begin to move when an infinitesimal extension 
of the crack releases more energy than is needed to create fracture surface. This idea 
is enormously successful in practice, but conceptually incomplete. We will show that one 
cannot understand how a crack moves without taking into account the fact that it moves in 
a fundamentally discrete and not continuous medium. In a lattice there are some velocities 
for which crack solutions do not exist at all, others where cracks are linearly unstable and 
accelerate to higher velocities, and others for which crack tips are unstable and break apart 
altogether. While these conclusions are all compatible with continuum mechanics, they 
were not predicted by it, and are somewhat surprising. 

The experimental observations which motivate this study are 
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• that cracks in amorphous brittle materials such as glass or Plexiglas pass al- 
most instantaneously from quasi-static motion to motion at about 10% of the 
Rayleigh wave speed (DOLL AND WEIDMANN 1976; TAKAHASHI, MATSUSHIGE, 
AND SAKURADA 1984) . 

• that they seldom travel faster than 60% of the Rayleigh wave speed (KOBAYASHI, 
OHTANI, AND SATO, 1974; RAVI-CHANDAR AND KNAUSS, 1984), although accord- 
ing to continuum theory (FREUND, 1990) this wave speed should be the limiting 

velocity. 

• that at about 40% of the Rayleigh wave speed, the acceleration of cracks 
slows sharply, their velocity begins to oscillate (FINEBERG et. al., 1991, 1992; 
WASHABAUGH AND KNAUSS, 1993) , they emit high-frequency acoustic waves 
(GROSS et. al., 1993) , the energy needed to form new crack surface increases 
dramatically (GREEN AND PRATT, 1974) , and the fracture surface shows periodic 
structure (MECHOLSKY, 1985) correlated with velocity osciUations (FINEBERG, 
et. al., 1992) . The basic time scale of the velocity oscillations is much larger 
than the time scale on which the crack snaps atomic bonds, and has been 
unexplained. 

The continuum theory of fracture contains many hints concerning these phenomena, but its 
predictions have always seemed ambiguous. In the very first detailed calculation concerning 
dynamic fracture, Yoffe (1951) showed that at around 60% of the Rayleigh wave speed, a 
crack should become unstable, since the maximum tensile stress would no longer be directly 
ahead of the crack, but would instead be off at an angle. The dynamical implications of 
this calculation were however unclear. Would the crack branch, would it begin oscillatory 
motion, or do something else? To make matters more puzzling, depending on precisely 
which component of the stress one monitors, the instability can set in at different velocities, 
and there is no criterion deciding between the various possibilities (FREUND 1990). For 
example, Freund (1974) has found that the velocity at which a moving crack consumes 
enough energy so that two nearly static cracks could travel instead of one fast one is 45% 
of the Rayleigh wave speed cr, while some recent perturbative calculations (XU AND KEER, 
1992; GAO, 1993) find a dynamical instability at around 65% of cr. In addition there are hints 
of instabilities in calculations coupling thermal and mechanical motion (LANGER, 1993A; 
LANGER AND NAKANISHI, 1993) , and in other simple energy balance arguments (SLEPYAN, 
1992; GAO, 1993; SLEPYAN, 1993) . 

All of these calculations are constrained by the fact that continuum elasticity is not 
well suited to describe the microscopic processes by which elastic energy is converted to 
broken bonds. We believe that understanding how this happens, even in the simplest case, 
is the key to resolving the qualitative puzzles mentioned above (MARDER AND LIU, 1993; 
LIU, 1993). 

The idealized models of brittle lattice fracture we will solve in this paper are not 
intended to describe material microstructurcs realistically, but, since they can be solved 
analytically, provide solid ground on which to show how the dynamics of fracture alter 
when discreteness is taken into account. In these models, steady state crack motion is 
impossible for velocities less than around 30% of the relevant wave speed, and at around 
50% of this wave speed, steady state crack motion becomes unstable. We will describe the 
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precise dynamical way the instability occurs, show some of its consequences after onset, 
and compare the results with experiment. 

The type of instability found in the lattice theory is similar to that known by the term 
intermittency (MANNEVILLE, 1990), and it occurs in a simple way. The crack tip starts 
from some initial condition, and tries to settle into a steady state configuration in which 
it snaps one by one the atomic bonds perpendicular to its direction of motion. Just as 
it is about to reach this state, a seemingly irrelevant bond parallel to the direction of 
motion snaps and throws the crack off its course. The crack tip is thrown far away from 
the steady state configuration, and gathers itself up to try the approach again, repeating 
the process periodically. One of the surprises in the analysis is that the time scale which 
governs the approach of the crack tip to the steady state is given by the time for transient 
perturbations to die out, and is formally just the inverse of the amount of dissipation in 
the model. Physical systems with small amounts of dissipation should be expected to show 
oscillations on relatively long time scales. This idea provides a tentative solution to the 
most perplexing experimental observation. 

From a formal point of view the calculations described here build upon the work of 
Slepyan (1981) and co-workers (KULAMEKHTOVA, SARAIKHIN, AND SLEPYAN, 1984) , lat- 
tice solutions for dislocations, (ATKINSON AND CABRERA, 1965; CELLI AND FLYTZANIS, 
1970; THOMSON, HSIEH, AND RANA, 1971) and upon results of prior numerical simulations 
(ASHURST AND HOOVER, 1976; SIERADZKI it et. al., 1988; MACHOVA, 1992) . An interesting 
new result is that many of the steady state configurations which have been derived previ- 
ously are unphysical and do not exist. In addition, the stability of the remaining steady 
states is examined carefully. The general rule is that any state whose velocity increases 
when one pulls on it harder is linearly stable; however, there are nonlinear instabilities to 
watch for as well, and the points where these occur are determined through a combination 
of analytical and numerical techniques. The smallest negative eigenvalue governing ap- 
proach of transient configurations to the steady state is calculated and found to be simply 
the inverse of the coefficient of dissipation in the model. 

There are some clear differences between the behavior of the lattice models and ex- 
perimental systems. The velocities at which various processes occur in the models are all 
distinctly higher than their supposed counterparts in experiment. Whether these differ- 
ences can be attributed to the fact that the model is schematic, or whether additional 
physical processes altogether actually operate in the experiment, will have to be deter- 
mined in the future. Some of the most important additional processes to consider have to 
do with the fact that the experiments are conducted in polymeric solids, not simple lattices, 
and are fully three-dimensional (RICE, BEN-ZION, AND KIM, 1994; PERRIN AND RICE, 1994) 
while the theory considers only two dimensions. The solutions considered in this paper do 
not allow dislocations (ZHOU, et. al., 1994) . This fact may curiously make comparison with 
polymers more appropriate than with crystals, despite the fact that the calculations are 
carried out in lattices. 

Because analytical calculations of crack motion in lattices are elaborate, the paper 
will proceed in steps. First, in Section 2 we will review the theory of crack motion in 
a continuum strip. Next, in Section 3 we will present a one-dimensional model which 
displays almost all the features of the more elaborate cases. In Section 4 we will solve 
a two-dimensional model in which mass points move with only one degree of freedom 
(Mode III). In Section 5 the mass points will be allowed to move with two degrees of 
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Figure 2.1: All the calculations in this paper will concern a semi-infinite 
crack moving through the center of an infinite strip. 

freedom (Mode I), and the the calculations compared with some new experimental results 
in Section 6. 



2. Crack motion in a continuum strip 

The successes of dynamic fracture mechanics (freuND, 1990) provide the backdrop to 
this paper. The most complete understanding of crack dynamics in a continuum has been 
obtained for a semi-infinite crack in an infinite plate, with the crack driven by loading 
on its faces (WILLIS, 1990) . However, we will study a different geometry here, that of a 
semi-infinite crack in an infinite strip (KNAUSS, 1966; BERGKVIST 1973; BERGKVIST 1974) 
. There are two reasons for the choice. First, the choice of a strip makes it simpler to 
compare with numerical simulations over long periods of time, since the effects of top and 
bottom boundaries are included in the calculation. Second, the strip geometry is the only 
one in which time-independent loading naturally leads a crack to move in steady state at 
a constant velocity. In an infinite plate, generic time-independent loading causes a crack 
to accelerate indefinitely. 

What is known about crack motion in a strip (Fig. 2.1)? In certain respects, the 
geometry is extremely convenient, since when a crack moves some distance I along the 
middle of a strip, all features of the problem are unchanged except that stressed material 
in advance of the crack has been converted to unstressed material behind it. For this 
reason, if elastic energy W is stored per unit length ahead of the crack, then the energy 
fiowing to a steadily moving crack tip to create new surface is also W per unit length. The 
expected asymptotic solution in a strip is that the crack will move at the velocity v such 
that 

w = r{v), (2.1) 

where r(f ) is the fracture energy per unit length crack advance. There is an apparent 
problem with this picture. What happens if T{v) is constant, but W is slightly more than 
F? In this event, the crack absorbs the excess energy by a continual slow acceleration. 
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Figure 2.2: Solution of Eq. (2.2) with F given by Eq. (2.3) and shown in 
inset, for Wq = 1, p = 1, cr = .9, w = 1, and v{t = 0) = 1 x 10^^. A long 
period in which the crack struggles to move is followed by a period of rapid 
acceleration, and a final approach to steady motion. 



moving adiabatically up through steady states and asymptotically approaching a terminal 
velocity given by a the Rayleigh wave speed (MARDER, 1991; LIU AND MARDER, 1991) . An 
approximate equation of motion which describes such cases is (LIU AND MARDER, 1991) 

v=i[l-T{v)/W]{l-v'/clf (2.2) 

with w the half- width of the strip, cj^ the Rayleigh wave speed, and c; the longitudinal wave 
speed. The right hand side of Eq. (2.2) is a simple approximation of elaborate analytical 
expressions; it is only accurate to around 10%, and that only when the dimensionless group 
vw/cf is small, but is very useful for qualitative analysis. 

The particular form of r(t') is not specified by continuum elasticity, and by varying T{v) 
one has the potential to find a wide range of fracture dynamics. Experiments show a very 
rapid jump in velocity at the onset of fracture, a fact we will explain later as a consequence 
of the lattice theory, but one can also explain it in the context of the continuum by studying 
a case in which fracture energy initially decreases with velocity, and then increases again. 
Taking 

rw=r(0)-M^„(^-(-^)^) (2.3) 

leads to the fracture dynamics shown in Fig. 2.2. States at f = for slightly less than 
r(0) are actually metastable, since crack motion would be possible if it could be initiated. 
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Figure 2.3: Solution of Eq. (2.5) witli r(0)/Wo = 1 and cr = 1. Only the 
solutions indicated by the thick line are stable. 



Another way to look at the case in which fracture energy T initially decreases as a 
function of v is to ask what steady velocities are possible as a function of the energy stored 
in the strip far ahead of the crack. In order to make a connection with results to be 
obtained in lattices, define 

A = ^W/T (0), (2.4) 

so that A is proportional to the vertical displacement of the sides of the strip ahead of the 
crack. Then solving Eq. (2.1) and Eq. (2.3) for v/cr gives 

V Wo ^ ' 16 4 Cfl' ^ ' 

which is pictured in Fig. 2.3. The point of this calculation is to show that if fracture energy 
is initially a decreasing function of velocity, there must be a forbidden band of steady state 
velocities; cracks which start on the lower branch shown in Fig. 2.3 quickly accelerate to 
the upper one, meaning that such fractures are activated. 

Thus within continuum theory one may have a range of velocities in which cracks are 
unstable and accelerate rapidly towards some minimum stable velocity. However, lattice 
theories predict something even more severe: that there is a range of velocities in which 
steadily moving solutions do not exist at all. 
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Figure 3.1: This one-dimensional model mimics the motion of a crack in 
a strip, incorporating effects of discreteness. One can view it as a model 
for the atoms lying just along the surface of a crack. The mass points are 
only allowed to move vertically, and are tied to their neighbors with springs 
which break when they exceed a certain extension. The lower portion of the 
figure shows an actual solution of the model, using Eq. (3.28), at d = 0.5, 
b=0.01. 

3. One-dimensional model 
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3.1. Definition and Energetics 

In order to mimic the motion of a crack in a strip, including the effects of underlying 
discreteness, but otherwise making the calculation as simple as possible, consider the model 
shown in Fig. 3.1. One can view it as a model for the atoms lying just along a crack surface. 
They are tied to nearest neighbors by elastic springs, with spring constant K — 1, and 
tied to a line of atoms on the other side of the crack line by similar springs, which however 
snap when extended past some breaking point. The lines of atoms are being pulled apart 
by weak springs of spring constant K = 1/N. These weak springs are meant schematically 
to represent N vertical rows of atoms pulling in series, and in later sections will be treated 
more realistically. The equation which describes the upper row of mass points in this model 
is 

Ufn+i,+ — ^Um,+ + Um-i,+ Elastic coupling to neighbors 

+jj {Um — Um,+) Driven by displacing edges of strip 

+ {um- - Wm,+) ^ (2w/ - \um- " «m,+ |) Bonds which snap 
—bum + Dissipation. 

(3.1) 

There are a few terms that need discussion. First, ^ is a step function, and the term 
containing it describes bonds which snap when their total extension reaches a distance 
2uf, where is a fracture distance. Second, a small amount of dissipation has been 
added, the term proportional to b. Originally this term was included for formal reasons 
in order to make Fourier transforms well defined, but it will eventually turn out to have 
physical importance. The amount of dissipation will usually be taken to be vanishingly 
small. Third, the height that mass points reach after the crack has passed is Un, and this 
term describes the force driving crack motion. 

In the discussion that follows, some of the equations will have boxes to the left, like 
Eq. (3.2). These are equations that are true regardless of the particular lattice that is 
being considered, and will be taken over without change in later sections. 

As varies, the natural scale on which a displacement Un is able to drive crack motion 
varies, so one should find the natural dimensionless constant which governs crack motion. 
The question to ask is, when is enough energy stored in the strip, per lattice spacing far 
to the right of the crack, to break one bond along the crack line? An important physical 
quantity to define in answering this question is obtained by going far to the right of the 
crack, and taking the ratio of the displacement of the atom just above the crack line, ?7right 
to the total displacement at the top of the strip Un. Denote this ratio by 

■ Qo ^ (3.2) 

Un 

Suppose that mass points far to the left and far to the right of the crack are stationary, and 
that dissipation is negligible. Far to the right of the crack, one has that ?7right = '^m+ = 
—Um- , and so balancing forces on masses with m » gives 

1 

2C/right = {Un - f^right) (3.3) 
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Remembering a factor of two for the upper and lower halves of the strip, the energy per 
bond for m is therefore 

4^ " ^"g^*)' + \ (2^"ght)' (3.5) 

^ £^right = 2Qo {Un? . (3.6) 
Far to the left of the crack, one simply has that 

■ C^ieft = Un. (3.7) 

The energy per bond for m <^ is just that needed to stretch the spring between Um+ and 
Um- from rest to breaking, and is just 

-E^break = 2uj (3.8) 

Assuming there is no other sink of energy, one must have that 

2Qo{UNf = 2u}. (3.9) 

Therefore, the minimTim value of Un at which enough energy is stored to the right of the 
crack so as to be able to snap the bonds along the crack line is 

■ Ufj = (3.10) 



For this reason, it is convenient to define a dimensionless measure of how far one has pulled 
the edges of the strip, 

Uf 

By definition, energy balance requires that steady crack motion is only possible for A > 1, 
and a system in which a crack moved when A = 1 would have to be perfectly efficient in 
turning potential energy into crack motion. 



3.2. Wiener-Hopf solution for steady states 

While the preceding calculations motivate the definition Eq. (3.10), they are based 
upon false assumptions. When a crack moves in steady state, Slepyan (1981) first showed 
that the mass points far from its tip are necessarily in motion. As a result the energy 
accounting carried out above is wrong, and crack motion is only possible for A > 1. The 
goal now is to examine steady states in detail. 

A steady state in a lattice is more complicated than one in a continuum; it is a configu- 
ration which repeats itself after a time interval 1/v, but moved over by one lattice spacing. 
In steady state, one has the symmetries 



(3.12) 
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Um+ (t) = U0+ (t - m/v) , (3.13) 

which means that all spatial behavior is contained in the time history of any single mass 
point. Take this mass point to be tto+, and denote it simply by u{t). Using Eq. (3.13) and 
Eq. (3.12) in Eq. (3.1) gives 

1 

N 



il = u {t - 1/v) - 2u + u {t + 1/v) + — {Un -u)- 2ud {2uf - 2\u\) - bu. (3.14) 



Eq. (3.14) can be solved analytically using Wiener- Hopf methods (NOBLE, 1959) . Here 
is the solution. 

There must be some first time at which u{t) rises to 1 and the 9 function vanishes. 
Let us take this time to be t = 0. Assuming that u rises above the height of 1 once and 
remains above it for good, one can write 

u = u{t- 1/v) -2u + u{t + 1/v) + ^ (^t/ive""l*l -uj - 2ue {-t) - hit. (3.15) 

Eq. (3.15) does not quite follow from Eq. (3.14), since the factor exp(— Q;|t|) has appeared. 
It is introduced just to make Fourier integrals converge, and a will tend to zero at the end 
of the calculation. 

Everything in Eq. (3.15) can be Fourier transformed in a straightforward way except 
for the term 6{t)u{t). Simply define 



u- 

and 



{u) = J dte''^^e{-t)u{t) (3.16a) 



U+ (cj) = J dte^'^^e (t) u {€) (3.166) 

so that U(u3)^ the Fourier transform of u{t), is 

U {u) = U+ (ou) + U- (a;) (3.17) 

The crucial observation is that is free of poles in the upper half complex uj plane, while 
U~ is free of poles in the lower half plane, since the integrals Eq. (3.16) are obviously 
convergent in these cases. 

Using these definitions, one can now transform Eq. (3.15) to read 

^/ / ...r . Un f 1 . 11 U 



with 



-uj^U = 2 (cosa;/-f; - 1) t/ + <^ + >-^- 2C/" + ioubU. (3.18) 

^U{oj)F (uj) - 2U- (a;) = -% I + ] (3.19) 

F (a;) = + 2 (coso^/w - 1) - ]^ + (3-20) 



Solving for XJ with the aid of Eq. (3.17) gives 



^+/,.^ , [1,1 
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Define 



F{u)-2 



and use the fact that a is vanishingly smaU, so one only needs the value of F{Q) 
on the right hand side of Eq. (3.21) to write 



1 1 

+ 



Q {u) U+ (uj) + U- (uj) = QoUn 

' a + lu) a — iuj 
Here Qq is given by Eq. (3.4), and one can easily check that 

■ Qo = Q (0) . 



(3.22) 

-1/N 

(3.23) 
(3.24) 



From a formal point of view, Eq. (3.23) is important because all of the lattice models 
in subsequent sections can be put in exactly this form, with the function Q becoming 
progressively more complicated as the model becomes more realistic, but everything else 
remaining precisely the same. 

The Wiener-Hopf technique directs one to write 



Q{u) 



Q H 

Q+ {co) ' 



(3.25) 



where is free of poles and zeroes in the lower complex co plane and Q"*" is free of poles 
and zeroes in the upper complex plane. One can carry out this decomposition with the 
explicit formula 



exp 



du' 
2^ 



e-*'^'Mng(w') 



(3.26a) 



exp 



lim 



duj' In Q ioj') 
27r iuj ^ e — iuj' 



(3.266) 



Now separate Eq. (3.23) into two pieces, one of which has poles only in the lower half 
plane, and one of which has poles only in the upper half plane: 

U+ (a;) 



QoUn 



1 



QoUn 



Q+(a;) Q- (0) {-ioj + a) Q- (0) {iuj + a) Q~ {u)' 



(3.27) 



Because the right and left hand sides of this equation have poles in opposite sections of 
the complex plane, they must separately equal a constant, C. The constant must vanish, 
or U~ and will behave as a delta function near t = 0. So 



U- {u) = Un 



QoQ (^) 
Q- (0) {a + iu) ' 



and 



U+ {io) = Un— 



Q- (0) {a - iu) 



(3.28a) 



(3.286) 
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One now has an explicit solution for [/(w). Numerical evaluation of Eq. (3.26a), and 
U{t) from Eq. (3.28) is fairly straightforward, using fast Fourier transforms. However, 
in carrying out the numerical transforms, it is important to analyze the behavior of the 
functions for large values of a;. In cases where functions to be transformed decay as 
this behavior is best subtracted off before the numerical transform is performed, with the 
appropriate step function added back analytically afterwards. Conversely, in cases where 
functions to be transformed have a step function discontinuity, it is best to subtract off the 
appropriate multiple of e~*6{t) before the transform, adding on the appropriate multiple 
of 1/(1 — ioj) afterwards. A solution of Eq. (3.28) constructed in this manner appears in 
Fig. 3.1. 



3.3. Relation between A and v 



There is an important point which has been forgotten. This solution is only correct if 
in fact 

u{t) = Uf &tt = (3.29) 

because this is supposed to be the moment at which the bond between uo+ and uq- breaks. 
The only parameter in the problem one is free to vary is Un, so the condition Eq. (3.29) 
chooses a value of Un, or its dimensionless counterpart, A. Once one assumes that the 
crack moves in steady state at velocity v, there is a unique A which makes it possible. 
To obtain Eq. (3.29), one needs to require that 

^e--^t/- (a;) = Uf. (3.30) 



lim 



This integral can be evaluated by inspection. One knows that for positive t > 0, 



I 



duj exp [— ■io't] U [uj) = 



(3.31) 



, and that any function whose behavior for large u is 1/iu! has a step function discontinuity 
at the origin. Therefore, Eq. (3.30) and Eq. (3.28a) become 



Uf = UnQo 



Q- (oo) 



(3.32) 



Q- (0) ■ 

Since from Eq. (3.22) follows that Q{oo) = 1, one sees from Eq. (3.266) that 

■ (oo) = Q+ (oo) = 1. (3.33) 

As a result, one has from Eq. (3.32) and the definition of A given in Eq. (3.11) that 

Q- (0) 



A = 



(3.34) 



To make this result more explicit, use Eq. (3.266) and the fact that Q{—lij) = Q{ui) to 
write 

' r dJl [lnQ(a;') , lnQ(-a;')" 
J 2^2 



Q- (0) = exp 



+ 



lUJ' 



e + ioj' 



(3.35) 



13 



exp 



2^ 



1 



-2iu' 



In 



Q{co') 



+ 



e2 + uj" 



Ing(O) 



■ ^Q- (0) = v^exp 

Placing Eq. (3.37) into Eq. (3.34) gives 

A = exp 



dio' 1 
27r 2iuj' 



In 



Q_(cy) 



(3.36) 



(3.37) 



(ia;' 1 
27r 2za;' 



In 



Q(a;0 



(3.38) 



This expression is not completely general because the fracture condition that U{t — 0) 
must equal Uf is not completely general. In the lattice considered in Section 5, one has 
instead that U{t = 0) must equal 2uf/y/3. Apart from this constant of proportionality 
everything goes through as above, and one has the general result that 



C exp 



27r 2iuj 



^|lnQ(a;') -lnQ(u;')} 



(3.39) 



where C is a constant of order unity that is determined by the geometry of the lattice, 
equaling 1 for the model of this section and the model of Section 4, and 2/y/S for the model 
of Section 5. When written in this form, Eq. (3.39) is suitable for numerical evaluation, 
since there is no uncertainty relating to the phase of the logarithm. 

When b becomes sufficiently small, Q is real for real uj except in the small neighborhood 
of isolated roots and poles that sit near the real uj axis. Let r^' be the roots of Q with 
negative imaginary part (since they belong with Q'^), r~ the roots of Q with positive 
imaginary part, and similarly pf^ the poles of Q. Then one can rewrite Eq. (3.39) as 



A = CJ JJ^^, for 6 = 



(3.40) 



Together with Eq. (3.28), Eq. (3.39) and Eq. (3.40) constitute the formal solution of 
the model. Since Q is a function of the steady state velocity v, Eq. (3.38) relates the 
external driving force on the system. A, to the velocity of the crack v. 

It is interesting to plot the function A{v) obtained from Eq. (3.40) (Fig. 3.2). Because 
all steady states occur for A > 1, one necessarily concludes that not all energy stored to 
the right of the crack tip ends up devoted to snapping bonds. The fate of the remaining 
energy depends upon the amount of dissipation b, and the distance from the crack tip one 
inspects. In the limit of vanishing dissipation 6, traveling waves leave the crack tip and 
carry energy off in its wake; the amount of energy they contain becomes independent of 
b. Such a state is depicted in Fig. 3.1, which shows a solution of Eq. (3.28) for v — 0.5, 
N = 9, and b = 0.01. For all nonzero 6, these traveling waves will eventually decay, and 
the extra energy will have been absorbed by dissipation, but the value of b determines 
whether one views the process as microscopic or macroscopic. 
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Figure 3.2: The velocity of a craclc f is plotted as a function of the 
driving force A for the one-dimensional model. The calculation is carried 
out using Eq. (3.40) for = 100 and v = 0.5. The thick upper line 
indicates physically realizable solutions, and the line along v — indicates 
the range of lattice-trapped solutions. 



3.4. Forbidden velocities 

The jagged structure of Fig. 3.2 makes it appear that many different states, emitting 
different quantities of radiation, can coexist for some values of A. This conclusion is largely 
incorrect, for two reasons. In Appendix III, it is shown that states are linearly unstable 
whenever v is a decreasing function of A. So all the backward-traveling portions of the 
curve can be ruled out. In addition a final condition has been neglected. Not only must 
the bond between wo+ and uq- reach length at t = 0, but this must be the first time 
at which that bond stretches to a length greater than Uf. For < f < 0.3 . . . (the precise 
value of the upper limit varies with b and N) that condition is violated. The states have 
the unphysical character shown in Fig. 3.3. Masses rise above height Uf for t less than 
0, the bond connecting them to the lower line of masses remaining however intact, and 
then they descend, whereupon the bond snaps. Since the solutions of Eq. (3.15) is unique, 
but does not in this case solve Eq. (3.14), no solutions of Eq. (3.14) exist at all at these 
velocities. There is a forbidden band of velocities. 

Nevertheless, multiple solutions for some values of A are still possible. The phe- 
nomenon of lattice trapping (THOMSON, 1986) allows a crack to sit still in a lattice under 
some range of external strains, before the first bond holding it snaps and it begins to 
extend rapidly. The lattice trapped solutions of this model are constructed in Appendix I, 
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Figure 3.3: The height of wo+ is plotted as a function of time for v — .3, 
b — .01, and = 9. Notice that u rises above Uf = 1 before t = 0, and is 
actuaUy descending at the moment when it again passes height Uf and is 
supposed to have broken. 



and shown to exist in the range 

Vs-i _ 



5176... < A < 



V3 + 1 
V2 



1.931 



These bounds do not correspond to the value of A obtained from Eq. (3.40) as v 
limit is carried out in Appendix II, and shown to be A = {y/b + l)/2. 



(3.41) 
0; that 



3.5. Linear Stability 

The stability of steady states can be studied in a straightforward manner, by adding a 
small term to them and linearizing in the perturbation. One finds that even for stable 
states, transients converge slowly to the final state, at rate e~^*, where b is the damping in 
the model. 

This final result can be established in a simple way which relies only upon time reversal 
symmetries as follows: 

Consider any equation of motion for some variables Umit) of the form 



Um (t) = Om {u) - bUm {t) , 



(3.42) 



where O is invariant under time translation, and even under time reversal. If one starts 
with a base solution u^{t — m/v), then perturbations of the form e^*'U^(t) = e^*u^{t — m/v) 
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will obey the eigenvalue equation 

q'^ul^ + 2qu]^ + = ■ - bu]^ (3.43) 

One easily checks that vP {t — m / v) is an eigenfunction with eigenvalue q — Q. In addition 
notice that 

u]^ = vP {-t - m/v) (3.44) 



obeys 



1 _ -A 



t) li 



■ + hiil^ (3.45a) 



do 

-2bul + ul = -^.u^- bill,. (3.456) 

Comparing with Eq. (3.43), one sees that to first order in b, one has an eigenfunction with 
eigenvalue —b, given simply by time reversing vP. 

This general conclusion is reproduced by the much more detailed analysis of Appendix 
IV, as discussed in Section 4.3. 



4. Simple Two-dimensional Model (Mode III) 



4.1. Definition and Energetics of the Model 

The calculation of the previous section will now be extended to a two-dimensional 
lattice model, depicted in Fig. 4.1. A crack moves in a lattice strip composed of 2(A'" + 1) 
rows of mass points. All of the bonds between lattice points are brittle-elastic, behaving 
as perfect linear springs until the instant they snap, from which point on they exert no 
force. The location of each mass point is described by a single spatial coordinate w(m, n), 
which can be interpreted as the height of mass point (m, n) into or out of the page. The 
index m takes integer values, while n takes values of the form 1/2, 3/2, . . .N + 1/2. The 
model is described by the equation 

u (m, n) — —bit + - J-' [u (m', n) — u (m, n)] , (4.1a) 

nearest , / 
neighbors ' 

with 

(u) ^ ue {2uf - \u\) (4.16) 

representing the brittle nature of the springs, 9 the step function, and b the coefficient of 
a small dissipative term. The boundary condition which drives the motion of the crack is 



u{m,±[N + 1/2]) = ±Un 



(4.1c) 
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0,N+l/2 




0, -N-1/2 



Figure 4.1: Lattice model of fracture. Tlie equilibrium locations of 
mass points are indicated by the white dots, while the black dots indi- 
cate the displacements u{m, n) of mass points out of the page once stress 
is applied. The top line of dots is displaced out of the page by amount 
u(m, N +1/2) = A\/2N + 1, and the bottom line into the page by amount 
u{m,—N — 1/2) = —A\/2N + 1. Lines connecting mass points indicate 
whether the displacement between them has exceeded the critical value of 
2 (see Eq. (4.16)) . The crack tip has just reached location m = 0. 



As before, it is important to find the value of C/jv for which there is just enough energy 
stored per length to the right of the crack to snap the pair of bonds connected to each 
lattice site on the crack line. For m 3> one has 

u{m,n) = — (4.2) 

so that the energy stored per length in the 2A'" + 1 rows of bonds is 
1 



X [2 Upper Bonds/Site] x [Rows] x [Spring Constant] xCZ/jgij^. 



2 



= -2(2Ar + l)- ^ (4.3) 

2 ^ ^ 2\N + I/2J ^ ^ 

= 2QQiUNf, (4.4) 

with Qq = 1/{2N + 1) given as before by Eq. (3.4). The energy required to snap two bonds 
each time the crack advances by unit length is 
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Therefore, as before, the proper measure of external driving is 

A UnVQo 



, (4.6) 

Uf 

a quantity which reaches 1 as soon as there is enough energy to the right of the crack to 
snap the bonds along the crack line. 



4.2. Reduction to Form of Previous Section 



The techniques used to solve this model were found some time ago by Slepyan (SLEPYAN, 
1981). There are differences between details of his solution and ours because Eq. (4.1) de- 
scribes motion in a strip rather than an infinite plate, and in a triangular rather than a 
square lattice. The strip is preferable to the infinite plate when it comes time to compare 
with numerical simulations, while reducing to the simpler infinite plate results in certain 
natural limits. 

Assume that a crack moves in this lattice in steady state, so that one by one, the bonds 
connecting u{m, 1/2) with u{m + 1, —1/2) or u{m, —1/2) break. They break because the 
distance between these points exceeds the limit set in Eq. (4.16) and as a consequence of 
the driving force described by Eq. (4.1c). Supposing that these are the only bonds which 
snap (an assumption to which we will return later) it is possible to calculate the motion 
of all the mass points above the line n = 1 /2 as a function of the mass points on the line 
n = 1/2, since in any region where the bonds do not snap the model has simple traveling 
wave solutions. 

In steady state, one has the symmetry 



u (m, n,t) — u{m + l,n,t + 1/v) 

and also 

u (m, n, t) — —u (m, — n, t — [1/2 — gf„] /f ) 
which implies in particular that 

u{m,l/2,t) = -u (m,-l/2,t- 1/2^;). 

We have defined 

' 



(4.7a) 
(4.76) 

(4.7c) 
(4.8) 



if n= 1/2,5/2... 
gn^{l ifn = 3/2,7/2..., 

mod (n — 1 /2, 2) in general 

One can now eliminate the variable m entirely from the equation of motion, by defining 



Un (t) = u (0, n,t) , 
and write the equations of motion in steady state as 

■+Wn+1 {t - {gn+1 - 1) /v) + Un+1 {t - Qn+l/v) 
-\-Un (t +l/v) - 6Un (t) +Unit- 1/v) 
,+Un-l {t - {gn-1 - 1) /v) + Un-1 {t - Qn-l/v) 



(4.9) 



1 



Un {t) = 2 

-bUn 



(4.10a) 
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if n > 1/2, and 

+%/2 (^) + «3/2 - 1/^^) 
+W1/2 (t + 1/v) - 4wi/2 (t) + Uy2 {t - 1/v) 
. + [u-l/2 {t) - «l/2 {t)] i-t) + [«_i/2 {t - 1/v) - Uy2 {t)] (1/ (2^) - t) _ 



1 

2 



(4.106) 



bu 



1/2 



The time at which the bond between tt(0, 1/2, t) and tt(0, —1/2, t) breaks has been chosen 
to be t = 0, so that by symmetry the time the bond between u{0, 1/2, t) and u{l, —1/2, t) 
breaks is 1/2^;. 

For n > 1/2 it is easy to solve the linear set of equations Eq. (4.10a). Fourier trans- 
forming in time gives 



-Un+i (a;) 
' 1 



^iu!{gn+i-l)/v _|_ ^ioj{gn+i)/v 



-u;\nH =ibuj+ +-Un{u;) e*'^/^ - 6 + e"*'"/^ 

^ L 



(4.11) 



Let 



(u) = uy2 (^) e'=("-V2)-i^ffn/(2^^). (4.12) 
Substituting this expression into Eq. (4.11), and noticing that Qn + Qn+i = 1 gives 

-■^^1/2 (^) Qi^{gn+l+gn-'2)/{2v) _|_ giw(5„+i+5„)/(2t;) 



^wi/2 (w) = «^>wwi/2 (c^) + -Wi/2 (w) M'^/'' - 6 + e-*'^/ 



+ -^«l/2 (^) e"'^ giw(s'„-l+fl'n-2)/(2t;) _|_ giw(c?„_i+5„)/(2v) 



cj^ + z6a; + 2 cosh [k) cos (a;/ (2v)) + cos (a;/!?) — 3 = 0. 



( 



.13) 



Defining 



z = 



one has equivalently that 
with 



3 — cos (w/v) — up" — ibio 
2 cos {uj/2v) 

y ^ z + \J z^ - 



y = e\ 



One can construct a solution which meets all the boundary conditions by writing 



Un (oo) = ui/2 {uj)e 



-y[N+l/2-n] _ y-[N+l/2-n] 



+ 



[/at (n- 1/2) 2a 



N 



(4.14) 

(4.15) 

(4.16) 
(4.17) 

(4.18) 
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This solution equals ^1/2 for n = 1/2, and equals UN2a/ {a^ + uP') for n = AT + 1/2, with 
a sent to zero at the end of the calculation. The most interesting variable is not wi/2, but 
the distance between the bonds which will actually snap. For this reason define 



U{t) = 



Ml/2 {t) - U-i/2 (t) _ Ui/2 (t) + Ui/2 {t + l/2v) 



Rewrite Eq. (4.106) as 



1 

2 



+'"1/2 + 1/^) - 4tti/2 {t) + Wl/2 {t - l/v) 

-2U (t) d (-t) -2U{t- l/2v) e{\/ {2v) - t) 



— hu 



1/2 



Fourier transforming this expression using Eq. (4.18) and defining 



(kue^^u (t) e (±t) 



now gives 



with 



u^/2 {u) F iu) - (1 + e^-/2-) U- (a;) 



Un 2a 
' N uj'^ + a'^' 



F{u) 



y[N-l] _ y-\N-l] 



yN _ y^N 



2z ) cos {ijj/2v) + 1 



Next, use Eq. (4.19) in the form 



(1 + 

U (a;) = ^ ^«i/2 (a;) 



to obtain 

Writing 
finally gives 

with 



Un 2a 
N uj'^ + a^ 



U (io) = U+ (io) + U- (io) 

U+ (uj) Q (uj) + U- (u) = UnQo + -^—] 

[a + iu; a — luj ) 

Q{io) = 



F{u) 



(4.19) 



(4.20) 



(4.21) 

(4.22) 
(4.23) 



(4.24) 

(4.25) 
(4.26) 
(4.27) 

(4.28) 



F{uj)-1- cos {ijo/2v) 

To obtain the right hand side of Eq. (4.27) one uses the facts that -F'(O) = —l/N, and that 
a is very small, so that the right hand side of Eq. (4.27) is a delta function. 
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Figure 4.2: Plot of U{t) for v ^ .5, N = 9, and b = 0.01, produced by 
direct evaluation of Eq. (3.28). Note that mass points are nearly motionless 
until just before the crack arrives, and that they oscillate afterwards for a 
time on the order of 1/b. 



Notice as promised that Eq. (4.27) is identical to Eq. (3.23). Therefore from this point 
forward the analysis of the Section 3 can be repeated without alteration. In particular 
Eq. (3.28) gives an explicit expression for U{(jj), and Eq. (3.38) and Eq. (3.40) describe the 
relation between A and v. 

Without repeating the calculations, we mention for later reference that a square lattice 
can be solved in the same manner - it is simpler than the triangular lattice - and the results 
are the same except that 



F{uj) 



y[N-l] _ y-[N-l] 



N 



replaces Eq. (4.23) and 



y 



Q = 



y 



-N 



-2z + l 



F-2 



(4.29) 



(4.30) 



replaces Eq. (4.28). 

There are three parameters to vary in looking for numerical solutions of Eq. (4.27). 
The most important is v, the velocity of the steady state. In addition, one can also control 
N, the width of the strip, and b the amount of dissipation. There is a natural limit in which 
many results become independent of N and 6, the limit of a macroscopic dissipationless 
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strip, namely N ^ oo and 6 — > 0. No physical results depend upon the order in which 
these limits are taken, although the integrals one has to perform for Eq. (3.38) look very 
different. For 1/b ^ N, InQ/Q is only nonzero near the finite number of points where 
Q has a pole or a zero. One can use Eq. (3.40) in this case. For N ^ 1/b these poles 
and roots merge into branch cuts, and the integrand of Eq. (3.38) becomes smooth. By 
taking the limit in this way, one can show the equivalence of results in the strip with 
the results found previously by Slepyan (SLEPYAN, 1981) for a square lattice occupying an 
infinite plate. In the following discussion, the calculations will be carried out for N = 9. 
Although seemingly far from the limit A/" — > oo , all physical quantities that have been 
checked so far change only in the fourth decimal place when N increases from 9 to 100. So 
= 9 is large enough to give an accurate picture of the behavior of the model, but small 
enough to make all types of computations rapid. 

All the phenomena discussed in relation to the one-dimensional model occur here. A 
picture of u{t) for v = .5 appears in Fig. 4.2, a graph of A(t') appears in Fig. 4.3, and 
a picture of an unphysical state at v = .2 appears in Fig. 4.4. For = 9 and b = 0.05, 
starting at the wave speed c = \/3/2 and working downwards, u{0) first becomes negative 
at V = 0.244, or 28% of the wave speed c. At lower values of v, there is no evidence that 
the steady states are ever physical. For example, at v = 0.106, ^(O) is positive. However, 
examining the state, one finds that u earlier rose above w(0), descended below it, and is 
on the rise again. Stationary states with v = are physical, and exist up until the point 
where they are met by the states with v 0. Once again, there is a velocity gap, and no 

steady states exist between v = and v = 0.244 Above v = 0.577, or v/c = .666, the 

steady state solutions are unstable, for reasons that will be discussed in a later section. 



4.3. Linear Stability 

The complete calculation of linear stability is carried out for the triangular lattice of 
this section and the details are relegated to Appendix IV. The results are as follows: 



Define 



2 

3 — cos {u>/v) + {q — iuj) — ibu 



(4.31a) 



2 cos {uj/2v) 




(4.316) 



and 




2zq > COS {u>/2v) + 1. 



it + 1/2^;) 



(4.33a) 



(4.32) 



and restrict attention to modes with the symmetry 




= —u 



-1/2 {t + 1/2^) , 



(4.336) 
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0.0 

1.0 1.5 2.0 

A 



Figure 4.3: Steady state velocities as a function of external strain A, 
calculated for strip half- width = 9 and dissipation & = using Eq. (3.40). 
The thick lines indicate cases in which the steady states are stable. Zero 
velocity states at strains A > 1 correspond to the phenomenon of lattice 
trapping (THOMSON, 1986). 



Then defining 



one has 

C/(a;) = (Q+-Q,^)^, (4.35) 

with 'uP{0) the velocity of the steady state u^{t) at t = 0, and Uq subject to the consistency 
condition 

Uo = J^U {u) , (4.36) 

which may be rewritten as 

1^(1- Q- (a;)) ^S, = u' (0) = So. (4.37) 

The spectrum of perturbations about steady states is given by the zeroes of Sq — Sq. If 
Sq — Sq has a root when the real part of q is positive, the corresponding steady state is 
unstable. 

One consequence of Eq. (4.37) is that any state whose velocity decreases as A increases 
must be unstable. This result may be established by considering the behavior of Sq for 
small q. If the slope of Sq is positive at the origin, then there must be a root for positive 
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5.0. 



4.0. 



3.0. 




Figure 4.4: Behavior of u{t) for v = .2, b = 0.05, and = 9. Notice 
that at t = 0, indicated by the dashed hne, u is decreasing, and that it had 
already reached height 1 ear her. This state is not physical. 



g, since Sq > 0, and S goes to zero as q 
expression, evaluated for the steady state: 



oo. The slope Sq is given by the following 



d 



or using Eq. (3.28), 



d\nA 
dliiv 



(4.38) 



(4.39) 



When accurate numerical evaluation of S{q) is necessary, it can be carried out from 
the integral 



Qq{co') 



(4.40) 



The moderately lengthy derivations of Eq. (4.39) and Eq. (4.40) are given in Appendix 

V. 

Fig. 4.5 shows a graph of Sq. The somewhat surprising triangular shape has a clear 
physical origin. According to Eq. (4.31a), when q is real and very small, it acts so as to 
shift the effective value of the damping b. Right at g = —b/2, the system is fooled into 
thinking that b has vanished, and as g decreases further, it is as if b has changed sign, or 
as if time has started to run backwards. The symmetry of Q implies that 



Q{u;,b) = Q{-u;,-b) 



(4.41) 
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Figure 4.5: The integral on the right side of Eq. {IV. IQ) is plotted as a 
function of g for d = 0.5 and h = 0.0125. 



Q- {uj, h) Q+ {-u, -b) = Q- {-uj, -h) Q+ {uj, b) 



(4.42) 



Since the left side of this equation has roots only above the real axis, and the right side 
only below, 

(4.43) 
(4.44) 



Q- (0,-6) = l/Q+ (0,6) 



so using 



(0) = 1 = V2iV + 1 



Q+(0)' 



one has that for small positive e 



V+e = -5^1ng+(0) = ^. 



(4.45) 



In short, the sign of the slope of Sq changes sign right at g = —6/2, resulting in a stable 
eigenmode with eigenvalue exactly —b. This result is exactly that obtained in Section 3.5. 
A picture of the mode is given in Fig. 4.6. 
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Figure 4.6: The eigenmode U{t) for q = b, given by u {—t) is pictured for 
V = 0.5 and b = 0.0025. 

4.4. Nonlinear Instabilities 

Showing that steady states are hnearly stable is not sufficient to demonstrate that they 
are physical. It was assumed in their derivation that the only bonds which break are those 
which lie on the crack path. From the numerical solutions of Eq. (3.28), one can test this 
assumption; it fails above a critical velocity. Recall that the sound speed c equals \/3/2. 
For A?" = 9, at a velocity of Vc/c = .666 . . ., Ac = 1.158 . . ., the bond between tt(0, 1/2) 
and tt(l, 1/2) reaches a distance of 2uf some short time after the bond between tt(0, 1/2) 
and w(0, — 1/2) snaps. The steady state solutions strained with larger values of A are 
inconsistent; only dynamical solutions more complicated than steady states, involving the 
breaking of bonds off the crack path, are possible. To investigate these states, one must 
return to Eq. (4.1) and numerically solve the model directly. These simulations have been 
carried out by Liu (MARDER AND LIU, 1993; LIU, 1993) , and some results are contained in 
Fig. 4.7. 

The diagram shows patterns of broken bonds left behind the crack tip. Just above the 
threshold at which horizontal bonds begin to break, one expects the distance between these 
extra broken bonds to diverge. The reason is that breaking a horizontal bond takes energy 
from the crack and slows it below the critical value. The crack then tries once more to 
reach the steady state, and only in the last stages of the approach does another horizontal 
bond snap, beginning the process again. This scenario for instability is similar to that 
known as intermittency in the general framework of nonlinear dynamics (MANNEVILLE, 
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A=1.835 
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Figure 4.7: Pictures of broken bonds left beiiind tiie craclc tip at four 
different values of A, from simulations of Liu (.) The top figure shows the 
simple pattern of bonds broken by a steady-state crack. At a value of A 
slightly above the critical one where horizontal bonds occasionally snap, 
the pattern is periodic. All velocities are measured relative to the sound 
speed c = Notice that the average velocity can decrease relative 

to the steady state, although the external strain has increased. As the 
strain A increases further, other periodic states can be found, and finally 
states with complicated spatial structure. The simulations are carried out 
in a strip with half- width iV = 9, of length 200 and b = 0.01. The front 
and back ends of the strip have short energy-absorbing regions to damp 
traveling waves. 

1990) ; the system spends most of its time trying to reach a fixed point which the motion 
of a control parameter has caused to disappear. 

Here is a rough estimate of the distance between broken horizontal bonds. Let Uh{t) 
be the length of an endangered horizontal bond as a function of time. Actually, one needs 
to view matters in a reference frame moving with the crack tip, so every time interval 1/v, 
one shifts attention to a bond one lattice spacing to the right. When A is only slightly 
greater than Ac, the length of such a bond viewed in a moving frame should behave before 
it snaps, as 

Uh -^2uc + ^{A- Ac) - Sue-^K (4.46) 

Here duh/dA means that one should calculate the rate at which the steady-state length 
of Ufi would change with A if this bond were not allowed to snap, and Su describes how 
much smaller than its steady-state value the bond ends up after the snapping event occurs. 
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Figure 5.1: Diagram of lattice used in this section 

From this expression, one can estimate the time between snapping events by setting to 
2uf and solving Eq. (4.46) for t. The result is that the frequency v with which horizontal 
bonds snap should scale above the critical strain Ag as 

a result that is consistent with the numerics, but hard to check conclusively. One can 
calculate numerically that duh/dA = 5.5 for the conditions of Fig. 4.7, but Su is hard to 
find independently. Assuming that Eq. (4.47) is correct, one finds from the second picture 
of Fig. 4.7 that 5u = 0.04. Further increasing the external strain A makes a wide variety of 
complicated behavior possible, including dendritic patterns, in the lowest panel of Fig. 4.7 
that we will see are reminiscent of experiment. 



5. Solution of Fully Two-Dimensional Model (Mode I) 

The triangular lattice solved in Section 4 contains all the physical phenomena we now 
know to be associated with lattice fracture. Still, it is not very realistic because of the 
restriction that mass points move in only one dimension. That restriction is relaxed in this 
section. Unfortunately, the calculations rapidly become so lengthy that they are difficult 
to explain. For this reason, most of the details are contained in a MAXIMA script which 
was used to produce the results. However, one can at least study the question of energy 
balance without too much trouble. 

The lattice studied in this section is depicted in Fig. 5.1. Each point has two degrees 
of freedom associated with motion in the plane. In the absence of external forces, the 
mass points all sit on a triangular lattice, and they interact linearly with their nearest 
neighbors. The interaction between two neighbors is a general linear function of both the 
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parallel distance between them, and the perpendicular distance between them. The two 
spring constants k_\_ and /cy are introduced so as to accommodate a general Poisson ratio. 
Choosing an arbitrary and k± results in long-wavelength transverse and longitudinal 
sound speeds 

3a2 



3a^ 
8m 



(3A;_L + , 



(5.1a) 



(5.16) 



with m the mass of the points, and a the lattice spacing. 

Unfortunately, for purely technical reasons, it is not possible to adopt this form of the 
nearest neighbor force everywhere on the lattice. In order to solve this model as in the 
previous sections, it is necessary that mass points on opposite sides of the crack experience 
a force only along the line between their centers. As partial compensation, however, one 
can take this spring constant k^^ to be different from the parallel force constant elsewhere 

in the lattice. Numerical calculations show that varying within reasonable bounds has 
very little effect upon the physical results. 

Taking bond interactions between neighbors in this way, one can form a correspondence 
with Section 3. Far ahead of the crack tip, the ratio of vertical displacements of points 
right above the crack line to vertical displacement at the top of the strip is 



1/4 (/cx + 3A;||) 



1/4 (A;^ + 3/c||) +2Ar (3/4) kf' 



(5.2) 



As a fracture criterion, one demands that bonds snap when the distance between neighbors 
across the crack line is 2uf greater than it is in equilibrium, so that the energy needed for 
fracture is 

^breat = S/cf^i (5.3) 



The displacement of the top of the strip just necessary to supply enough energy to snap 
these bonds is given again by Eq. (3.10) and Eq. (3.11). 

The main difference between the calculation here and the one in Section 3 results from 
a consideration of the fracture criterion. Letting •u(0, 1/2) be a point on the crack line, its 
bond with a neighbor on the other side will snap when 



\u{0, 1/2) -u{0, -1/2) I > 2uf 



(5.4) 



^ (0, 1/2) - uy (0, -1/2)) + i {u- (0, 1/2) - (0, -1/2)) 



>Uf. 



(5.5) 



However, if one searches for the variable in terms of which Eq. (3.23) and its sequels 
remains true, one finds instead 



{uy (0, 1/2) - uy (0, -1/2)) + („- (0, 1/2) - u^ (0, -1/2)) 



(5.6) 
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Figure 5.2: Velocity v measured in units of the transverse wave speed ct 
as a function of A, for the model of this section, calculated for = 15, 6 = 
0.012, kj_ = —.3, = 2.666, = 2.366. The thick lines are physically 
realizable states, the thin solid lines are linearly unstable, and the dotted 
lines are unphysical states. 



is what is needed. The reason is that far ahead of the crack, U defined this way approaches 
■u^, and all the displacements far to the right are purely vertical. In terms of this variable, 
fracture occurs when 



U{t = 0)=uf 



(5.7) 



as claimed before Eq. (3.39). A second difference has to do with the form of the dissipation. 
Experimentally, waves of frequency uj decay at a rate 6(a;) ~ oj. This does not seem well 
understood theoretically, but in any event, the calculation uses dissipation in the form 



(5.8) 



rather than iiohu as before. The constant is chosen to be something smaller than 
frequencies of interest. 

What remains is to calculate the function Q appropriate to this lattice. It is a lengthy 
task, carried out with symbolic algebra, and relegated to Appendix VI. No really new ideas 
enter; some of the results are summarized in Fig. 5.2 and Table 6.1. Fig. 5.2 indicates that 
there are small separate bands of low-velocity states, separated by small unstable regions. 
Note, November 1995: Mode I fracture is subject to several instabilities beyond those studied 
here. The complete story has yet to be worked out. 
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6. Comparison with Experiment 

From the lattice theory emerge three predictions which can be tested experimentally. 

1. For a band of velocities beginning at zero and proceeding up to around 30% of the 
relevant wave speed, steady crack motion is impossible. 

2. Steadily moving cracks emit radiation. 

3. At around 60% of the wave speed crack tips become unstable, repeatedly attempting 
to branch, and creating microcracks whose spacing is governed by the attenuation rate 
of high-frequency sound waves. The critical velocity for this instability depends on 
details of microstructure. 

We have carried out experiments in Plexiglas (PMMA), differing from those previously 
reported (FINEBERG, et. al, 1991, 1992) because of improvements in instrumentation that 
lead to a five-fold increase in the accuracy of velocity measurements. The experiments 
are conducted in thin strips whose geometry is chosen to match as closely as possible con- 
ditions of the theory. Results from the various lattice models are compared with these 
experiments in Table 6.1. All of the phenomena seen in the lattice are present in the 
experiment, although the velocities at which they occur are different. Since changing 
from a triangular to a square lattice changes the critical velocities substantially, this dis- 
crepancy with experiment is not surprising. Plexiglas is a polymeric solid, and bears no 
microstructural resemblance to a triangular lattice. It is certainly possible to study more 
complicated lattice theories numerically in search of quantitative agreement, but in the 
following discussion, we want mainly to emphasize the qualitative correspondence between 
the experiments and the theory. 

The first prediction has been verified whenever crack dynamics have been measured 
carefully in brittle materials (TAKAHASHI, MATSUSHIGE, AND SAKURADA 1984) , although 
not generally given much significance. In PMMA, the jump from quasi-static to rapid 
motion goes to a velocity of around 175m/s, which is 18% of the Rayleigh wave speed; 
three of our experimental runs at three levels of externally imposed strain are shown in 
Fig. 6.1. The large accelerations seen in this figure are not inconsistent with continuum 
theory, as one sees by comparing with Fig. 2.2. From Eq. (2.2) with q = 2000m/sec 
and w = A cm, continuum theory predicts accelerations on the order 10^ m/sec^, or 10^ 
m/sec//i-sec. However, it is difficult within a continuum framework to account for the fact 
that cracks can be made to initiate for such a wide range of external conditions, and that 
the velocity at which rapid acceleration abruptly terminates is fairly independent of the 
amount of energy pouring into the crack tip. From the perspective of the lattice theories, 
both these facts are natural. Crack motion is an activated process, and can therefore begin 
for a range of external conditions, while the rapid acceleration ends when the crack velocity 
has passed through the forbidden band. 

Conceptually, one likes to think of cracks becoming unstable when their infinitesimal 
extension releases more energy than it consumes, so that the crack slowly accelerates up 
to high velocities. Unfortunately these slow extensions are dynamically forbidden. The 
crack must begin its dynamic life at high velocities, and the criterion for fracture initiation 
must be understood as a nucleation event. These remarks would seem to contradict the 
fact of quasi-static crack motion. In the context of the lattice models, the contradiction 
can be resolved if quasi-static growth is always a thermally or chemically aided process. In 
addition, apparent quasi-static motion may actually be a stick-slip process, in which a crack 
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Lattice 






One-d (AT = 9, 6 = 0.01) 


0.3 


- 


Mode III, triangular(7V = 9, 6 = 0.01) 


0.28 


0.67 


Mode III. square(A^ = 9.6 = 0.01) 


0..36 


0.76 


Mode I {N = 15) 


0.33 


0.56 


Mode I {N = 30) 


0.33 


0.55 


PMMA 


0.2 


0.33 



Table 6.1: Calculations of the minimum nonzero steady-state crack ve- 
locity, t'niirij and the velocity Vc at which the branching instability begins 
for several lattices. In all cases, the calculations are carried out in strips 
2A'" + 1 lattice points hight. Mode III is the lattice of Section 4, and Mode 
I is the lattice of Section 5. In the One-d and Mode III cases, velocities are 
measured as fractions of the sound speed, while in the Mode I case, and in 
PMMA they are measured as fractions of the transverse wave speed. The 
Mode I lattice parameters are chosen to fit the Poisson ratio, and measured 
ultrasonic attenuation of PMMA, and are k± = —.3, = 2.666, k^^ = 2.366; 
experimentally dissipation is observed to be proportional to frequency, and 
this scaling is used in the calculation as well, with the experimental coeffi- 
cient of 0.012, (JACKSON, PENTECOST, AND POWLES, 1972) . The compar- 
ison has been made with the plane strain Poisson ratio; to compare with 
the plane stress Poisson ratio, one can use = instead of k± = —.3, but 
the results are not substantially different. 

alternates between rapid motion and arrest on scales too fast for ordinary measurements to 
detect. This observation may have implications for fracture testing. Recent measurements 
of the fracture energy of PMMA have all been in the range of 200 to 350 J /m^ (kATSAMANIS 
AND DELIDES, 1988) . Our experiments involve long center edge cracks in long strips, so that 
the fracture energy can be deduced simply by measuring the stress per area in the strip. 
Inspection of Fig. 6.1 shows that the fracture energy of our samples can be as low as 80 
J/m^ - three times less than typical previous measurements. A^ote, November, 1995: These 
low fracture energies have not been reproduced in subsequent experiments, and this claim is 
incorrect. Fracture energies are consistent with previous measurements. One gets a sense 
of how close one is to the minimum fracture energy by watching the rate of acceleration 
after the crack initiates. Wc suspect that many of the experimental techniques used to 
measure fracture energies involve dynamical effects at a more important level than has 
been appreciated, but this point requires additional investigation. 

The second of the predictions has not really been verified experimentally. It is true 
that moving cracks always emit sound (GROSS, et. ai, 1993). However, acoustic transducers 
detect only up to around 10 MHz, which is far below the frequencies at which atomic 
bonds snap when a crack moves at hundreds or thousands of meters per second. At such 
high frequencies, one should expect sound to thermalize rapidly, and manifest itself as 
heat. Certainly there is always a large temperature rise near crack tips (GREEN AND 
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Figure 6.1: Three measurements of velocity versus time in Plexiglas 
(PMMA). The experiments are conducted in strips, and the energy stored 
per unit area to the right of the crack is printed each picture. Cracks are 
made to initiate with different energies available per length by preparing 
the notch which initiates fracture in different ways. Even if the system 
is stressed so gently that the crack does not accelerate noticeably after it 
begins rapid motion, there is still a jump in velocity up to around 200 
m/s. The upper part of the figure should be compared with Fig. 2.2. The 
fracture energy of the sample depicted in the lower part of the figure is 
several times lower than any other value for PMMA in the literature; see 
Katsamanis and Delides, (1988). 



PRATT, 1974; FULLER, FOX, AND FIELD, 1983) , but there are many potential sources for 
it, especially plastic deformation. To the extent that the high-frequency sound decays 
within the core region surrounding the crack tip, it causes no conceptual difficulty for the 
continuum picture of fracture. Our prediction is that the size of the heated zone should be 
set by the dissipation coefficient b, but this estimate has not been tested experimentally. 

The third of the predictions corresponds to observations in PMMA, which observe the 
emergence of microcracks after about 40% of the Rayleigh wave speed. The first publication 
we are aware of describing these as an important factor in the fracture of PMMA is by 
Doyle (l983); their role has also been emphasized by Ravi-Chandar and Knauss (1984). Our 
use of thin strip samples enables us to obtain cracks traveling in steady state at different 
mean velocities, and locate the onset of micro-cracking with greater precision. A drawing 
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of these cracks soon after onset is displayed in Fig. 6.2. We can use Eq. (4.47) to estimate 
the spacing of these branches. 

The main experimental ingredient needed for the estimate is acoustic attenuation at 
high frequencies, which has been measured in PMMA by Jackson, Pentecost, and Powles 
(1972). They find that attenuation per cycle is nearly constant in PMMA as a function of 
frequency and temperature, and obeys 

aA = 0.02 - 0.07, (6.1) 

where a is attenuation per length, and A is the wavelength of the waves one wants to 
consider (the experimental measurements are for longitudinal waves.) The lower value of 
0.02 holds up to 100 MHz; there is a gap in the measurements, and 0.07 is measured at 10 
GHz. Taking a\ — 0.02, one has that attenuation per time h is given by 6A/c/ — 0.02, with 
Q the longitudinal sound speed. Polymer units of size 0.5 [x separate coherently in PMMA 
so PMMA will be treated as a lattice of masses separated by a =0.5/um. Now the waves 
being excited by the passage of the crack have wavelength A = aq/f , where a =0.5//m is 
the lattice spacing and v is the crack velocity, so 6 ~ 0.02v/a ~ 12 MHz, for a ~ 0.5//m 
and V ~ 300 m/sec. Taking [9'U/i/9A](A — Ac)/5u/i to be on the order of 10%, one has 
from Eq. (4.47) and Eq. (6.1) the crude estimate that microcracks appear with a frequency 
of 

6MHz=^ Spacing =0.05mm (6.2) 

for a crack traveling at 300 m/sec. Close examination of the 1/32 inch samples of PMMA 
reveals a comb of micro-cracks, shown in Fig. 6.2, resembling those in Fig. 4.7, and ex- 
tending into the surface with a typical separation of about 0.07 mm. This spacing is larger 
than the rough theoretical estimate, but of the same order of magnitude. Both in the 
simulations and in the experiment, microcracks on the upper and lower portions of the 
surface are staggered. 

There is also experimental evidence for the claim that the onset of the microcracking 
instability is sensitive to details of lattice structure. As shown in Table 6.1, moving from a 
triangular to a square lattice pushes onset of the instability to higher velocities; the general 
lesson should be that cleavage along weak lattice planes discourages micro-cracking. In- 
deed, cracks reaching the transverse wave speed (and in one case claimed to exceed it) have 
been measured in cleavage of LiFl (OILMAN, KNUDSEN, AND WALSH, 1958), while CottereU 
(1964) has measured speeds approaching the Rayleigh wave speed for cracks traveling along 
prescribed grooves in PMMA. 

The comparison between simulation and experiment has to be made with two qual- 
ifications. First, the simulations are on strips only nine atoms high, from all points of 
view much smaller than the experiments. Second, the experiments are clearly not two- 
dimensional. The cracks seen in Fig. 6.2 are those that were visible within a particular 
plane, about 100 n thick. Changing the vertical plane changes the details of the pattern 
of microcracks; this observation partly explains why the microcracks in the figure appear 
intermittent. 

In connection with this three-dimensional structure, wc note that Perrin and Rice 
(1994) have recently found an interesting implication of elastic theory. When a crack front 
is considered in three dimensions, it is linearly stable against local fluctuations in fracture 
toughness, but upon repeated contact with inhomogeneities, the crack front can deform 
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Figure 6.2: Drawing of small microcracks that begin to emerge from a 
crack after it passes a velocity of around 330 m/sec in Plexiglas (PMMA). 
The basic geometry is the same as that shown in Fig. 4.7, with a spac- 
ing on the order of 0.07 mm. This diagram shows the microcracks near 
the point where they first appear; later, they are longer and more densely 
packed. The microcracks shown in the figure are those that appear in a 
microscope whose depth of field is about 100 fi, and the pattern of microc- 
racks changes as one looks at different vertical planes by changing the focus 
of the microscope. 



considerably, and experience velocity fiuctuations on the order of a quarter of the Rayleigh 
wave speed. Therefore, a crack moving nominally at one-third the wave speed may locally 
be moving much faster (or much slower). This picture provides one possible explanation 
for the fact that microcracks appear experimentally at velocities noticeably lower than 
that predicted by the lattice models. 

In prior publications (fineberG et. al, 1991, 1992) , our group has described a 500 KHz 
frequency emitted from PMMA, corresponding to a 1 mm scale structure on the surface. 
This spatial scale is several times larger than the one we have attributed to the onset of 
microcracking; the 1 mm oscillations begin in earnest at a somewhat higher velocity (400 
m/sec) than the initial microcracks (330 m/sec), and correspond to bunches of especially 
dense microcracks grouped at roughly 1 mm intervals. We do not yet know whether this 
phenomenon is the natural evolution of the microcracking instability when driven to large 
amplitudes, or whether new physical processes are coming into play. Thus the original 
oscillatory instability with which we began our investigations is still not well understood. 

While resolving numerous questions that arise in continuum models, and providing a 
qualitative understanding for many features of the experiments, the lattice models raise as 
many questions as they answer. When should one get micro-cracks, and when dislocations? 
How do dynamic ductile-brittle transitions occur? What are the effects of thermal noise? 
How is quasi-static motion possible? What would happen in a random lattice? How should 
one treat a polymer? How can the results be generalized to three dimensions? How does 
the instability progress towards macroscopic branching? What happens in larger-scale 
simulations with realistic bond forces? What are the implications for fracture testing? 
These are some of the points that deserver further investigation. 
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Appendix I. 

Calculation of Lattice- Trapped States 



In this Appendix, we will calculate the static solutions of Eq. (3.1). Dropping the 
subscript +, one begins with the equation 

= Um+l - 2Um + Um-l + ^ {Un - ^m) - '^UmOm, (-^-1) 

where 

/, / 1 for m > / r o^ 

L for m < 

For this solution to follow from Eq. (3.1), one must have that 

uq < 1 and U-i > 1. (-^-3) 

In other words, the bond at zero has not snapped, so it must be stretched less than a 
distance of 1, but the bond at —1 has snapped, so it must be stretched more than 1. The 
solution of Eq. (/.I) is 

_ ( Ur (to) = 2N+1 ^" '^rUl!^ for TO > 



\ ui (to) ^Un+ uiyf' for TO < ' ^ 
where 

yi-2 + 1/yi - 1/N = 0, (/.5a) 

and 

yr-2 + 1/yr - 2 - 1/N = 0. (/.56) 

In order for the left and right solutions to match onto one another smoothly, one must 
demand that 

ui (0) = Ur (0) and ui (-1) = Ur (-1) . (1.6) 

The four equations Eq. (/.5) and Eq. (/.6) determine ui, Ur, yi and yr as functions of Un. 
Using Uc = \/2N + 1, turning the condition Eq. (7.2) into a condition on A, and working 
in the limit N ^ oo gives 

A < }^I±1 for uo<l, {1. 7a) 

v2 

and 

A > y^zl for u-i > 1 {I.7b) 
v2 

These are the boundaries for the region of lattice trapping. 



Appendix II. 

Evaluation of A for small and for large v 



In this Appendix, we will evaluate A(w) for the one-dimensional model in the limits 
V — s> and — > 1, starting from Eq. (3.40) 
It is possible to evaluate the product 



where for generality Wi is the root of 

a;2-4sin2 {uj/2v) - AA^ 



(//.I) 



(//.2) 



analytically in certain limits. In the limit of low velocity, it is helpful first to write the 
condition that Eq. {II -2) have a root as 



uj = 2v sin ^ y^a;2/4 - A^. 



(11.3) 



In the limit of low velocity, by looking at a graph of Eq. {1 1. 2), one sees that the roots are 
given approximately by 



= 2v 



w- = 2v 



(no + z) TT — sin ^ ( \J [no + z]^ tt^ — A? 



(no + z) TT + sin ^ { \j [no + z]^ tt^ — A? 



(IIAa) 



(IIAb) 



and that there is always one more positive root than there are positive roots . The 
starting integer no is given by 

A = 2vno'K {II. ba) 

(one will have to fiddle around with v a bit to make this precisely true) and the largest 
value of z is nj, given by 

\/ A^ + 1 = (no + n/) VTT {IIM) 

(this comes close to the truth in the limit of small v). 

For the moment, let us restrict ourselves just to the positive roots of Eq. {1 1. 2). Then 



\\nP{A)^ ^In 



(no + ?') ^ ^ sin 



-1 



t'^ [no + 



Uf — l 

El" 



i=o (no + i) TT + sin ^ (^-sj v"^ [no + if tt^ — 

1 — sin"-^ v"^ [no + if tt^ — A'^^ / (no + 
1 + sin~^ ( \J [no + if tt^ — j / (no + i) tt 



+ In [(no -I- n/) tt] 

(//.6) 

+ ln[(no + n/)7r] (//.7) 
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Uf — l 



(—2) sin ^ I [no + if tt^ — A'^j / (no + i)7r + ln [(no + nf) n] 

J— n V / 

-2\ sin~^ Vx"^ - 



/ dx 

J A 



and integrating by parts gives 



X 



1 12/" 
-lnP(yl) = ln- + - / 

2 w TT 



X In xdx 



which after the change of variables x"^ = A^ + sit? 9 becomes 



P{A)^^ exp 



- / In [^2 ^ sin2 6*1 
ttJo 



de. 



In the particular case of Eq. (3.40), one wants to evaluate 



exp 



1 r, reV4+ l/4 + sin2^^ 
In 



e2/4 + sin^6' 



Gradshteyn and Rizhyk (1980) have the integral (4.399) 

'1 + yi + a 



dO. 



J In (l + a sin^ x) = 2n In ^- 



This gives finally 



A = 



v/iV4TT74 + aAT?74TT74 



e/2 + Vl + e2/4 



In the limit e — > 0, one has that 



A = i±^ = 1.6180..., 



(//.8) 
(//.9) 



(I/.IO) 



(//.II) 



(//.12) 



(//.13) 



(//.14) 



(//.15) 



the golden mean, in agreement to three places with the direct evaluation of the roots in 
Eq. (3.40). One has that a stationary lattice crack in a noiseless environment will not 
begin to move until the driving strain exceeds by this amount the strain that would be 
predicted in a continuum model. 

At velocities that approach 1, there is one root /+ and one root g'^. In the limit of 

smaU e, and for v = 1, g+ = 1.91892. One finds /+ = ^2^36 + 12(1 - v). Therefore 



1.91892 



\ ^2V3e + 12(l-^;) 



(//.16) 



Appendix III. 

Calculation of A for small damping 



This Appendix describes the calculation of A(v) in the limit of zero damping, 6 = 
for the triangular lattice of Section 3. Begin with Eq. (3.40). 

In order to make use of this expression, one must find an efficient way to locate the 
roots and poles of Q. This task is equivalent to finding the roots of F, defined in Eq. (4.23), 
and the roots of 

^(a;) = F(a;) - 1 -cos(a;/2v), {IH-'^) 

which is the denominator of Q. Both F and G have poles, but at the same values of u, 
and these poles cancel. The roots fall into two classes. In the first class arc roots for which 
1^1 > 1. These roots are well separated, and can be located by a standard root-finding 
algorithm. In the second class are roots for which \z\ < 1. These roots are very finely 
spaced, and for each region where \z\ < 1, there are of order N roots. These may be 
located by looking for the values of y where F vanishes. These are denoted by y^, and are 
indexed by j, which varies from 1 to AT. Inverting F to find y one has 



F 

Vj =e 



i(2j-l)7r/(2Ar+l)+ln 



Vj -cos(aj/2t)) 



/(2Ar+l) 



(m.2) 



Therefore 



F 

Zj = COS 



(2j-l)7r 1 f yf -cosu!/2v 

2N + 1 ^z(2iV+l) Uf (cosa;/2^;)-l 



(///.3a) 



Similarly, inverting G to find y and z one gets 



Zj = cos 



2j7r 
2N+1 



{IIIM) 



The roots of G are therefore very easy to find. By finding solutions of Eq. (Ill.Sa) one 
finds the roots of F. Once all the roots are located, one determines the sign of their 
imaginary part for infinitesimal b, and finally inserts them into Eq. (3.40) to find A. 



Appendix IV. 

Linear stability of steady state solutions 



This Appendix carries out the hnear stabihty analysis of steady state solutions for the 
Mode III model in detail. To begin, it is useful to recall the symmetries of the steady state 
solutions. These are 

u^{m,n,t) ^u^ {m + l,n,t + l/v) {IV.la) 



and also 



u 



(m, n, t) = -u^ (m, -n, t - [1/2 - g-n] /v) 



(IV.lb) 



Therefore the perturbations can be taken first to be eigenfunctions of the operator which 
translates m by 1, and t simultaneously by 1/v, perturbations with time dependence of 
the form e'^*f{n,t — m/v). Second, the perturbations are eigenfunctions of the operator 
which inverts states around the x axis while moving t forward by l/2v. If this operation is 
repeated twice, it simply produces one of the translations of the first operator, and must 
therefore have eigenvalue e^*. The operation carried out only once must therefore have 
eigenvalue ±e^*/^; even and odd modes, associated with the inversion symmetry. 
Therefore the states will be of the form 



iIV.2) 



1 
2 



(IV.Sa) 



(n, t — m/v) + {n, t — m/v) e^*, 

with v} small, and with both even and odd inversion symmetries possible. 

Placing this state into Eq. (4.1a) and expanding to first order, one has that 

(f'u^ (m, n) + 2q'u} (m, n) + (m, n) = 

^-v} (m + ^n+i - l,n+ 1) +tt^ {m-\- gn+i,n-\- 1) 
-\-v} (m — 1, n) — Qv} (m, n) -\- v} {m-\- 1, n) 

{m + gn-i - l,n - 1) (m + gn-i,n - 1) 
—hv} (rn, n) 

if n > 1/2, taking Uf = 1 

(f'v} (m, n) + 2qv} (m, n) + v} (m, n) = 

(m + gn+i - l,n + 1) + {m + g-n+i, n + 1) 

(m — l,n) — 4u^ (rn, n) + {m + 1, n) 

+ [u^ (m, n-l)-v} (m, n)] [9 {-t) - 6 {l - (t))] 

_+ [u^ {m + l,n - 1) - (m, n)] [9 (l/2v - t) - S {l - (t - 1 /2v))] 
— hv} (m, n) 

if n = 1/2. The boundary condition is that vanish at n = ±(A^+ 1/2), and is defined 
by Eq. (4.19) as the diff'erence between two mass points above and below the crack line in 
the steady state. Upon Fourier transforming Eq. {IV.2>) one sees that it differs from the 



1 
2 



{IV.Zh) 
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solution of Eq. (4.1) only in two ways. First, —uP' — > (9 — %oS)^ ^ changing the definition of 
z in Eq. (4.15). Second, the driving force on the crack face is no longer the applied stress 
(Too, but is now a set of two delta functions appearing near the end of Eq. (/F.36). 
Fourier transforming Eq. (/K3) one has therefore that 



Zn = 



3 — cos (co/v) + {q — too) — ihu 
2 cos {00/ 2v) 



yq = Zq+ \ Z\ - 1, 



{IVAa) 



{IVAh) 



and 



- ^ 1 die-* { Y-1I2 it) - <I2 (^)] ^ (-^) + [«-i/2 - 1/^) - «!/2 (^)] ^ (1/2^ - ^)} 



«-l/2 (0) - «}/2 (0) 



+ e 



ioj/2v 



u\i^{-\l2v)-u\i^{\l2v 



1/2 



2w0 (0) 



with 



(/F.5) 
(/F.6) 

At this point, one needs to distinguish between the modes which are odd and those which 
are even under inversion. Define 



[iV-l] _ -[iV-l] 



U{t) = 



w}/2 W - «il/2 (*) 



1/2 



for the modes such that 
and 

for the modes such that 



u\/2 (t) = -«-i/2 {t-l/2v). 

y ^ ^V2 (0 - ^-1/2 (t) 
«}/2 (t) = «?_i/2 - 1/2^) , 



(/K7a) 

(/V.76) 

(/K8a) 
(IVM) 



The least stable mode is of the type U{t), which will therefore be of interest in what follows. 
One has that 



or 



U (a;) Fg-{1 + cosuj/2v) t/" = - (1 + cosa;/2v) 



{IV.9) 



(IV.IO) 



45 



with Qq given by Eq. (4.28); the subscripts q are meant to remind that z depends upon q. 
This may be rewritten 

Qt Qq \Qq Qq) «° W 

This expression can be decomposed immediately to give 



Qq{u) \Qq Qq{00)J (0) ' 

where a constant term has been included to avoid a delta function singularity in U (t) at 
t = 0. Since as before Qq{(jo) goes to one at when uj goes to infinity 

U- (C) = (1 - Qq) {IV.U) 

Because of time translation invariance, one expects there to be a zero mode, whose 
eigenfunction is vP{t). By comparison with Eq. (3.28), one sees immediately that this is the 
case, since Qq=o appearing in Eq. (/y.l3) is identical to the function Q defined previously 
by Eq. (4.28). The general consistency condition which tells which values of q are allowed 
is 

u-(t = <» = }^l (1 - Q,- (-)) (IV.U) 
^ (0) = / ^ (1 - Q,- M) ■ (iVM) 

Using the fact that that the zero mode is just u^, one can also write 

0=1^ (Qo (^) - Qq (^)) =Sg- So. {IV.16) 

When this equation has a solution for q with positive real part the corresponding steady 
state is unstable, and otherwise it is stable. 



Appendix V. 

Cracks which go faster when pulled harder are stable 



This Appendix derives relation Eq. (4.39). 
Notice from Eq. (IVAa) that 

Zq {auj, av, b/a, q) — Zq {u, v,b,q — iau) + iui) 

Therefore 

Qq (acj, af , 6/a, q) = Qq (a;, v,b,q — iau + iu) 
(since the additional dependence of Qq upon uj and v is of the form oj/v) 

dq du> iu) dv iu db 

Now 



Q" (cj) = exp 



J—oo J 



so 



dq 

which implies, since Qq depends upon q only through z, that 



Considering the terms in Eq. {V.5) one by one, we have first 



/_ 



—oo 




r dte'^^ f —e-"^'H{it)lnQq 
J-oo J 27r 



' — oo 

duj J_oo J 27r 



d 

^^InQg (a;) 



Second, 



/° dte^'^^ / ^e-^'^'^J^lnQ, 



lim [ -T——, :r- — :lnQ„ 
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n / 



lim , „ 

= li„ 



1 1 

i{uj — uj') iu' 



Putting together Eq. {V.Q) and Eq. (V.l) gives that 



(K7c) 
{V.ld) 



d_ 
dq 



— InQ^ {uj) = 



.8 V d ^ b d 
du iuj dv iu> db 



In 



Qq (^) 
(0) 



(K8) 



Returning with this result to Eq. (4.37) gives 

S'o- 



Since 

one has finally that 



5^ 



ainA 
91nf 



{V.9) 



(V.IO) 



(V.ll) 

(V.U) 



An additional expression makes it possible to evaluate Sq more quickly, and without 
needing recourse to Fourier transforms. One has that 



lim iuj < exp 

w— >oo 



dte 



tcut 



duj' 
2^ 



- 1 



(F.13) 



(K14) 



The large a; behavior in the exponential is given by the discontinuity in the function 
of t at 0, so one has next 



lim iu { exp 

w— >oo 



1 f duj 



ioj J 2tt 

du' 
2^ 



InQg (a;') 
InQ, (a;'). 



- 1 



One has finally that 



ZTT 



Qqiu') 



{V.15) 

(F.16) 
{V.17) 



Appendix VI. 

Calculation of Q for fully two-dimensional model 



This Appendix describes the calculation of Q for the model of Section 5. There seems 
no simple way to report it, and the calculations were all carried out in MAXIMA; the 
batch file is available upon request from the first author. 

Let Aj describe the displacements between points as shown in Fig. 5.1. The restoring 
force parallel to the direction of equilibrium bonds will be while that perpendicular 
to this direction will be k^. The force due to the displacement of the particle along 
Ai = Wj-ij+i 



Ui^j IS 



{VI.l) 



-1 /-I V3 



2 ' 2 



^1) 



Vs 1^ 

T' 2 



^Af , -A^ 
2 ^' 2 1 



{VI.2) 



Adding up contributions from other particles in this way we get for the force due to 
neighbors 



D 

F (m, n) = ^ ^ kqdqi [Ki (m, n) ■ dqi^ 



iVI.3) 



In steady state, for n > 1 /2 the forces become 



1, 3, 



<+i (t - i9n+i - 1) a/v) + <+i {t 
+<-i (t - (fn+i - 1) a/v) + <_i {t 



Qn+ia/v) 



u 



n+1 



Qn+ia/v) 

{t - {9n+i - 1) a/v) - ul^-^ (t - Qn+ia/v) 

1 {t - {9n+i - 1) a/v) + u\^^ {t - Qn+ia/v) 



- ^< (t) 



+ku [< {t + a/v) + < (t - a/v) - 2< {t)] 



1, 3, 

4^^ + 4^1 



K+i - {9n+i - 1) a/v) + ul^^ [t 
{t - {gn+i - 1) a/v) + {t 



py = 



Qn+ia/v) 

Qn+ia/v) 

Qn+ia/v) 

Qn+ia/v) 



(VIAa) 

- 4< (t) 



<+i {t - {9n+i - 1) a/v) - <+i {t 
{t - {gn+i - 1) a/v) + <_i {t 
+k± [ul {t + a/v) + ul{t- a/v) - 2< (t)] 

(VIAb) 

In Eq. (VIA), all the displacements u are evaluated at m = 0, and this index is therefore 
dropped, so that the only remaining index refers to the layer number n. For n = 1/2, one 
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has instead 



1, 3, 

4^11 + 4^^ 



<+i - {9n+i - 1) a/v) + <+i (t - Qn+ia/v) 
+e i-t) «_i (t - (gn+i - 1) a/t;) - < (t)) 
+^ - t) «_i (t - Qn+ia/v) - < (t)) 

<+i {t - {gn+i - 1) a/v) - w^^i (t - gn+ia/v) 



2< (0 



(-t) «_i (t - (^n+i - 1) a/v) - ul (t)) 
+6 {a/2v - t) «_i {t - gn+ia/v) - (t)) 
+kn [< {t + a/v) + < (t - a/v) - 2< (t)] 



1, 3, 

4^^ + 4^1 



(F/.5a) 
2< (t) 



Wn+l - - 1) ^/'") + ^1+1 - 9n+\a/v) 

+e i-t) «_i (t - {gn+1 - 1) a/^) - < (t)) 
(a/2^; - t) (t - gn+ia/v) - < (t)) 

= r <+i (t - - 1) a/v) - (t - ^n+ia/t;) 

+^ - ^) -0 (t) «_i {t - (^7n+i - 1) a/v) - < (t)) 
_ +e {a/2v - t) «_i {t - gn+ia/v) - < (t)) 
+k± K + 0/-^) + < - aZ-y) - 2< (t)] 

(VIM) 

Now take the Fourier transform in time of these equations. For the layers with n > 1/2 



one has 



py 



1, 3, 

4^11 + 4^^ 



+k\\ 



1, 3, 

4^^ + 4^1 



+k± 



-<-i 



,iw(p„+i-l)a/v _|_ gioipn+ia/wj 
^giw(g„+i-l)a/t; _|_ g(i-ff„+io/t;)^ 



Juj{gn+i-l)a/v _ iLjg, 



(VI.Qa) 



gZw(5„+i-l)a/-(; _|_ ^iu)g„+ia/vj 

^iLj{gn+i-l)a/v _|_ g(t-5„+ia/t;)^ 

^iu){gn+i-l)a/v _ ^iu)g„+ia/v'^ 



HAj{gn+i-l)a/v _ iujg, 



{VI.6b) 



V _j_ ^—iuia/v 2 



50 



Substituting in the form 



{VI.7) 



gives 



+ 3/cj_) cos {u>a/2v) 



{y + y-') 1 



- {moj'^ + iLoh (a;)) Ux = \ +2k\\ cos (oja/v) - 3 + 



(Vl.Sa) 



-\/3i {k± — /c||) sin (a;a/2u) 



{y-y ^) 



U, 



(/c_L + 3/c||)cos(a;a/2^;) 



(moj^ + iojb (oj)) Uy=y +2k_i cos (ua/v) -3{k_\_ + k 

—VSi (/c_L — sin {ua/2v) 



(VIM) 



{y-y ^) 



Ux 



Here, and in what foUows, the dissipation coefficient b is understood to be of the form 



b{uj) = bo^Lj'^+Lol. 
The condition that the determinant of this system vanish determines y by 



{VI.9) 



+ 3/cj_) cos {u)a/2'i 



{y + y~^) 



Q _ -\-muP' + iLob + 2/c|| cos {ua/v) — 3 + 
+ 3 



(/cx — sin (ioa/2v) — — ^ — ^ 



It is more convenient to define 



z = 



2/ + 1/2/ 



(A;i + 3A;||) cos (u;a/2t;) '-^^ ^ ^ 
+moj'^ + iLob + 2/cj_cos {ioa/v) — 3 + /cy) 

(F/.IO) 
(F/.ll) 



Then the determinental condition becomes 

A = 3 - V le/c^/cil cos {auj/ {2v)f 

B = cos {aixj/2v) (2 (^3A;| + 2k\\k±_ + 3A;^) cos (acj/v) + 4 + (ma;^ + ia;6 - 3 + 

C = (ma;^ + iub — (fcy + (3 — cos (acu/f )))^ — (/cy — /cj_)^ ^cos {ato/v)'^ + 3 sin (acu/ (2f ))^ j 



z± = 



-B ± - AAC 
2A 



(F/.12) 
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There are now 4 values of y which satisfy Eq. (F/.IO) for any given a;, namely 



y± = z± + y {z±) - 1 



{VI.13) 



and two others given by the inverse of these, or equivalently by changing the sign of the 
square root. 
Define 



ma;^ + icob + (/cii + 3A;j_) cos {uja/2v) 



{y^ + y±) 



+2A;|| cos {(jja/v) — 3 + 
E± = {k± - kn) sin {ioa/2v) ^ '^^ ^ 



(VI.U) 



Then a general solution of Eq. {VI. 6) is 



V 



-iojgn/'2v 



(n-1/2) 

y+ 



, (n-l/2) 

+y- 



+ Un 



[n - 1/2) fo 
N 



E. 



(-n+l/2) / E+ 



D 



+ 



(-n+l/2) / E 



U2+ 



U2- 



(F/.15) 



The four functions ui±,U2± can be determined from the four conditions 



(F/.16) 



and 



u 



1/2 

y 

1/2. 



E+ 
D+ 



U2+ + 



E- 



Ul- + 



E- 



U2- 



(V/.17) 



Once they are determined, one can obtain in particular a solution for ^3/2) in terms of it^^g 
and i^[i2^ however, the expressions are too long to list here explicitly. 

Thus the problem is reduced to that of finding u^i2 and w^^g- These are determined 
by taking the Fourier transform of the equations on the line n = 1/2. Unfortunately, we 
do not know how to solve the equations in full generality. The one restriction that must 
be imposed is that right on the crack line, we must take k± = 0. Otherwise the formalism 
dies. We would like very much to overcome this restriction, but do not now see how to 
do it. However, to try to make up for it, we will let kn equal some arbitrary kj, on the 
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interface. Given this restriction, Eq. {VI.6) becomes 

4/2 (t) + 4/2 (t - «/^) - 2</2 (t)] 



+ 



1, 3, 



i-t) (u^_,/, it) - (t) 

+A0 ^a/2v - t) (^^1/2 {t - a/v) - uf/2 (t)) 
-kje {a/2v - t) (uy_y^ {t - a/v) - vF^^ (t)) 

</2 + «/^) + </2 - «/^) - 2</2 (0 



1, 3, 

4^^ + 4^1 



^/2 (*) + ^/2 - ^1'") - 2^1/2 (^) 



3fcf 
4 



^ (-t) [u\i^ (t) - u\i^ it)) 



+ 



x/3 



- t) {u\i^ {t - a/v) - u\i^ (t)) 
- /cy) [4/2 (0 - 4/2 - «/^) 

-/c||^ (a/2^; - t) («^i/2 - a/v) - u\i^ it)) 
'"1/2 + fi/'J^) + ■J^i/2 {t — a/v) — 2'u-^y2 (^) 



(y/.18a) 



The important property of this set of equations is that there is really only one linear 
combination of and which multiplies the 9 functions. This combination is 

(0 = ^ [4/2 + «/2^) - 4/2 {t)] + \ {4/2 + «/2^) + 4/2 (^) } (^^-19) 

To see why it enters, notice that when the strip is loaded in Mode I, one must have 
the symmetries 

4i/2 {t) = -4/2 + «/2^) {VI. 20a) 



4i/2 (t) = 4/2 (t + «/2^) 



{VI.20b) 
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This symmetry allows one to eliminate the fields with subscript —1/2 from Eq. {VI.18). 
The result is 



-k 



«3/2 (^) + «3/2 - a/v) - 2u\,^ (t) 



''3/2 



1/2 



px _ 



\/3 r 



4 11/ 1^3/2 "3/2 

</2 + (^h) + </2 - «/^) - 2</2 (0 

[C/ {t) e {-t) -U{t- a/2v) e {a/2v - t)] 



8VS 



{VI.21a) 



py 



«3/2 (^) + ""3/2 - - 2«l/2 (*) 



\/3 r 

+— (/e^-Zeil) «^/2 W - 4/2 - «/^) 



"l/2 + ^Z"^) + '"'1/2 - a/v) - 2'U^/2 (^) 



-4 (t) ^ (-t) + U{t- a/2v) 9 {a/2v - t)] 
8 

It is now possible to Fourier transform Eq. (y/.21). The result is 

4/2 {u;) (l + e-«/" 



1, 3, 

4^11 + 4^^ 



+^{k^-k 



"3/2 ('^) ( 1 



^iwa/v 



+k\\ 



2u^^2 COS (auj/v) — 2uy2 {^) 

u- H I 



8^3 



^iuja/2v 



{VI.21b) 



2</2 (^) 



(V/.22a) 



(i^^ + i^ll) K2M (l + e^^«/'')-2«^/2(-^) 



^4 
4 



+ ^(/c^-/c||)[4/2M (l-e-«/«)^ 



2u^y2 {^) COS {ato/v) — 2u^^2 {^) 



k^ 

{VI.22b) 

Using Eq. {VI. 15) to find andttg^g; Eq. {VI. 19) to find t/ (a;), one can eliminate 
aU variables but U from Eq. (F/.22). 
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Analyzing the a; — > behavior, where 

U u\i^, (V/.23) 

^3/2 ^ (1 - V^) «i/2 + f^iv/iV5 (o;) 

one finds finally that 

Q (a;) t/+ + t/- = QoC^iv<^ (a;) (F/.24) 
as in Section 3. The expression for Q is again too lengthy to record here. 



